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Abstract 

It is shown that the inverse Lagrangian map for the solution of the Burg- 
ers equation (in the inviscid limit) with Brownian initial velocity presents a 
bifractality (phase transition) similar to that of the Devil's staircase for the 
standard triadic Cantor set. Both heuristic and rigorous derivations are given. 
It is explained why artifacts can easily mask this phenomenon in numerical 
simulations. 



1 Introduction 

There is a renewed interest among physicists and mathematicians in the Burgers 
equation 

d t v + vd x v = vd\v. (1) 

In particular it has been discovered that, when the initial velocity vq(x) is a Brownian 
(or fractional Brownian) motion function of the space coordinate x and the limit of 
vanishing viscosity v is taken, the Lagrangian map 

a I— > x(a, t) (2) 

from the initial fluid particle position a to its position x at time t is a Devil's staircase 
P], 0, Pf ; this has consequences for the distributions of masses of large-scale structures 
in the Universe 0. 

Since the Lagrangian map (Q) is nondecreasing (a consequence of the fact that 
fluid particles can merge but not cross), two nonnegative measures may be defined via 
their increments. The direct Lagrangian measure associates to a Lagrangian interval 
[a, b] the length Ax = x(b) — x(a). The inverse Lagrangian measure associates to an 
Eulerian interval [x,y] the length Aa = a(y) — a(x). When the initial density field is 



1 



uniform, this length Aa is proportional to the mass in the interval and will therefore 
be called the "mass" . The multifractal properties of the direct and inverse Lagrangian 
measures can be analyzed by studying the scaling properties of the moments of their 
increments. Here, we shall be interested only in the inverse Lagrangian measure, 
whose moments are defined by 

M q (l) = ([a(x + l)-a(x)] C! ), q>0, (3) 

which does not depend on x because the Lagrangian map has homogeneous (station- 
ary in the space variable) increments. 

In Ref. it was conjectured that the inverse Lagrangian measure has a bifractal- 
ity similar to that known for the standard Devil's staircase associated to the triadic 
Cantor set (see Fig. [1]). According to the conjecture, the scaling exponents r q , ob- 
tained from the small-/ behavior of the moments 

M 9 {1) oc V\ (4) 

should present a phase transition: for q > q+, one should have r q = 1, while for 
< q < q*, one should have r q = q/q*. Phase transitions of this sort are frequently 
observed in studying fractal sets of physical interest (see, e.g., Ref.^j). Preliminary 
numerical tests, reported in Refs. [3] were rather inconclusive. Further (unpub- 
lished) simulations indicated the presence of a small-/ scaling regime with r q = 1 for 
any q. While developing the correct theory for the inverse Lagrangian map, we found 
that this is actually an artifact inherent to discrete numerical simulations which can 
hide the true scaling. 

Here, we demonstrate the conjectured bifractality. Section |2| is devoted to recall- 
ing some known results about the solution of the Burgers equation and to formulating 
the problem with (fractional) Brownian initial conditions. In Section |3] we present a 
heuristic approach (the physicist's viewpoint); the arguments are not rigorous, but 
encompass both the Brownian and the fractional Brownian case; for pedagogical rea- 
sons, the arguments are presented first for the standard Devil's staircase. In Section f|, 
we present a rigorous proof for the Brownian case. In Section [5], we present numerical 
simulations, discuss the aforementioned artifact and show how to overcome it. 

2 The Lagrangian map for the Burgers equation 

We recall the construction of the solution to (0) (see Refs. [|, || for details). The 
solution at time t, for a continuous initial velocity Vq(x), is given by 

v(x,t) = v (a(x,t)). (5) 

Here, x i— > a(x, t) is the inverse Lagrangian map obtained as follows : we denote by 
^oO^) the initial velocity potential (v(x) = —d x ip(x)); then the Lagrangian location 
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a(x,t) is such that ipo{ a ) — {x — a) 2 /(2t) achieves its global maximum. This map is 
nondecreasing and has discontinuities at shock locations. Its inverse a i— *> x(a, t) is 
called the Lagrangian map. Those Lagrangian locations which are not mapped into 
shocks are called regular. 

For the case of fractional Brownian initial velocity, vo(x) is a random function, 
defined from — oo to +00, which is Gaussian and satisfies the following conditions 
(angular brackets denote averages) : 

(vq(x)) = 0, (6) 
([v (x') - v (x)} 2 ) = C 2 \x'-x\ 2h , (7) 

with < h < 1. The case h = 1/2 is the standard Brownian motion curve (in the 
space variable) for which it was shown in Ref. [l[] that the regular points form (almost 
surely) a set of Hausdorff dimension D = 1/2. For the case h 7^ 1/2, very strong 
numerical evidence was given in Ref. || that D = h. It follows that the Lagrangian 
map is a Devil's staircase. 

The scaling properties of the fractional Brownian motion Vq(x) and of the Burgers 
equation imply that |3|] 

v (x, t) '= t^v (xt~i=z, l) , (8) 

where = denotes probabilistic equality in law. Hence, the knowledge of the statistical 
properties at time t = 1 gives those at any other t > by simple rescaling. Henceforth, 
in discussing the theory, we shall sometimes set t = 1 without loss of generality. In 
performing numerical simulations with a finite mesh and a maximum box size L, the 
choice of the time t becomes relevant since the only scales which can be meaningfully 
related to the untruncated problem are between the mesh-size and L.f\ By changing 
t we can adjust the coalescence length 

l c (t) = (Ct) 1/{1 - h) (9) 

so that it does not lie too close to either of these. Over separations / ^> l c (t) the 
velocity increments, as well as the Lagrangian measures, remain basically unaffected, 
because not enough time has elapsed for substantial particle merging. In other words, 
the Lagrangian map is close to the identity. 



3 A heuristic approach 

We begin by recalling in some detail how bifractality arises for the standard inverse 
Devil's staircase (Fig. |I|). Let a denote the coordinate on the vertical (Lagrangian) 

1 Or a subset thereof, as we shall see in Section 0. 
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axis and x the coordinate on the horizontal (Eulerian) axis. At all (dyadic) x's which 
are integer multiples of 2~ n there is a jump (shock) of the inverse Lagrangian map 
a(x). It is straightforward to show that there are 

2°jumps of amplitude 3 _1 

2 1 jumps of amplitude 3~ 2 

(10) 

2 n jumps of amplitude 3~ (n+1 ) 

(This formula, which gives the distribution of the shock amplitudes, will be referred 
to as the mass distribution.) 

What is the equivalent of the quantity M q (l) defined above ? First, we shall restrict 
ourselves to increments / between two successive dyadic points x n>p = p2~ n and x njP+ \ 
of the same nth generation, so that I = 2~ n . The average (Q) of a quantity Q will be 
just a sum over the 2 n intervals of this type, divided by 2 n . For example, 

M q {l) = (Hx + l)-a{x)] cl ) (11) 

will be evaluated for I = 2~ n as 

p=2 n -l 

M 9 (J) = (l/2») £ [a{x n , p+l )-a{x n , p )f. (12) 

p=0 

Now, an important remark: If we do not restrict the Eulerian positions to dyadic 
points Xn tP , or if we work with a randomized version of the Devil's staircase, the 
Lagrangian increment a(x') — a(x) will be a sum over all the Lagrangian (shock) 
intervals corresponding to the dyadic Eulerian points between x and x' . This will 
generally include infinitely many intervals of generations n' > n. However, since the 
number of such intervals grows as 2 n , while their length decreases as 3~ n , the sum is 
dominated by just the first few terms for which n' is equal to n, or just a bit larger. 
Hence, we can replace the Lagrangian increment over an Eulerian distance / = 2~ n by 
just the contribution coming from the nth generation, thereby committing an error 
which just affects constants but not scaling. We thus have 

p=2 n -l 

M q (l) ~ (l/2 n ) £ [length of pth Lagrangian interval of nth generation] q . (13) 

The mass distribution gives us the number of shock intervals of a given length. 
Summing over all the generations m < n, we obtain : 

m=n 

M q {l) ~ (l/2 n ) ]T 2 m [3- (m+1) ] 9 . (14) 

m=0 
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It is clear that we must distinguish two cases : 

(i) when q > D = In 2/ In 3, the sum is dominated by its first term and we obtain 

M q {l) ~ (1/2") ~ l\ (15) 

so that r q — 1; 

(ii) when < q < D, the sum is dominated by its last term and we obtain 

M g (l) ~ (l/2 n )2 n [3- (n+1) ] 9 ~ 3- nq ~ Fln3/ln2 ; ^ 

so that r q = gin 3/ In 2, which establishes the bifractality for the standard Devil's 
staircase. This is seen as a rather elementary result, obtained without recourse to nu- 
merical computations, unlike the much more intricate phase transitions encountered 
in the study of certain circle maps ||| . 

We now turn to the Burgers problem for arbitrary < h < 1. The mass function 
is known [1], [| : the mean number per unit length of Lagrangian intervals with a 
length between Aa and Aa/2 is 

N(Aa) oc (Aa)- h . (17) 

Sinai's theory of the mass function |j] was done for t — 1. It is useful to introduce the 
correct dependence on the time of N(Aa,t). This is done by nondimensionalization : 
we multiply N(Aa, t) by the coalescence length l c ~ tVC 1 -^ and we divide the length 
Aa by l c . We thereby obtain : 

N(Aa,t) ~ {Aa)~ h /t. (18) 

For the standard Devil's staircase, the length of the Eulerian interval / was chosen 
as the inverse of the mean number (per unit length) of Lagrangian intervals for the nth 
generation. For the Burgers case, taking Aa = 2~ n and using (|18|) , this mean number 
is 2 nh /t. This is of course also the mean number of shocks per unit (Eulerian) length 
with a mass of the order of Aa = 2 _n . Hence, the corresponding Eulerian interval, 
the mean distance between two such successive shocks, is 

I ~t2~ nh = t(Aa) h . (19) 

This may also be interpreted as the length of an Eulerian interval such that the mean 
number of shocks having a Lagrangian length between 2~ n and 2 -(n+1 ) is order unity. 

We can now just repeat for the Burgers problem essentially the same argument 
as developed for the standard Devil's staircase, to obtain : 

i m=n i m=n 

M *(M) ~ — h E N(2- m ,t){2- m ]« ~ i £ 2 mh [2- m Y , (20) 

A m=0 1 m=0 
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Figure 1: Inverse Devil's staircase and the corresponding increment function 



where the upper limit n is given in terms of / by (|K|). 

As before, depending on whether q > q* = h or < q < q+, the sum is dominated 
by its first or its last term. Hence, we find, in the former case, that M q (l, t) oc I 1 and, 
in the latter case, that M g (l,t) oc l q l } \ so that r q = 1 for q > h and r q = q/h for 
< q < h,. This establishes the bifractality. 




4 Bifractality of shocks: rigorous theory 

Let a(x) be a (random) nondecreasing function defined on a finite interval J C R. 
Denote by M q (a(-), I) the space averaged value of the g-th power of the increment : 



M q (a(-)J):=—— [a(x + l) - a(x)] q dx } (21) 

-/ £ J I — 

where I — (I) means that the integration is over the interval I with the exception of 
an interval of length I at the right end of the interval J, where the increment is not 
well defined. 

Our aim is to study asymptotic properties as I — > of the mathematical expecta- 
tion (mean value) of M q (a(-),l). 

In Section |4.1| we shall derive estimates at first for the inverse of the standard 
Devil's staircase and then for a general function of this type, assuming that we know 
the statistical properties of N^, the number of shocks of magnitude of order 2~ k . 
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In Section fO] we prove inequalities for Nk for the case of the solution of Burgers 
equation with a Brownian initial velocity. 



4.1 Estimate of M q (a{-), I) 

Let a(x) be the inverse Devil's staircase function on the interval / = [0, 1] and let 
a n be its approximation constructed in the same way as for the n-th step of the 
construction of the standard Cantor set (see Fig. p. To calculate the functional 
M q (a(-), I) we calculate at first the sequence of approximations M q (a m (-), 2~ n ), m>n 
for the functions a m and then investigate its behavior as n — > oo. 

For any fixed n < m the increment function a m (x + 2~ n ) — a m (x) is a piecewise con- 
stant function with 2 k pieces of length 2~ n of magnitude 3~ k + 3 _n for k — 0, 1, . . . , n 
(see Fig. [I]). Therefore, 

l r 1 -^ 



M q (a m (-), 2- n ) := — / [a m {x + 2~ n ) - a m {x)] q dx 

1 — 2 Jo 

-i n— 1 

= y 2~ n 2 k (3~ {k+1) + 3- n Y 

1 - 2- n V 1 

k=0 

-i n— 1 i n— 1 

= y 2- n 2 k 3- ( - k+1)q + y 2- n 2 k 3- {k+l)q ■ R k , (22) 



k=0 k=0 



where 



Rk := (1 + 3-+ fe +Y - 1 < { g^T+r 1 ' || I < J! ■ (23) 



The first term in (|22|) can be calculated as follows : 

1 „,„,,^„ 2- n 3- g ^ / 2 \ fc 2-™3- g l-fi 



y 2- n 2 k 3~ {k+1)q = - — - — y [—) =— (24) 

while the second term becomes negligible compared to the first one as n — > oo. 



Consider now a more general function a : / — > R of this type. For each integer 
m > we construct an approximation a m (x), which is a piecewise linear function 
with Nk jumps of magnitude 6 [2~ fc_1 , 2~ fe ) for < k < m. Then the integral of 
the corresponding increment function with the increment 2~ n may be estimated from 
above by the integral over a piecewise constant function with 2 k pieces of length 2~" 
of mag nitude 2~ k + 2~ n for < k < n. There is a similar estimate from below (with 
the mag nitude 2- fc - 1 + 2" n ). Notice that our estimates do not depend on m, provided 
m > n. 



7 



Lemma 4.1 Let the function a m be as above and let us assume that C 1 2 fc//2 < Nk < 
C2 k ' 2 . Then 

2- w ! 4 ^ <MJa m n,2-)< °^ „ _ 2- 1 



IJI _ 2 -n 1 _ 2-(9~l/2) " 9V mw ' I- 1/1 _2-« 1 _ 2-(9-V2) ' 

Proof. In the same way as we did it for the inverse of the standard Devil's staircase 
we can estimate M q (a m (-), 2~ n ) as follows: 

M 9 (a m (-),2~ n ) := 1 / [a m (z + 2-") - a m (x)]« 



|/| - 2~ n i/-(2-«) 
2— n n— 1 



fc=0 

2~ n C n-1 



< Z 2 (fc+1)/2 (2~ (fc+1) I 2~ n ) g 



fc=0 



/m q— no — (g— 1/2) n— 1 

E2-^ 1/2) (l + ^), (25) 



I -2"" . , 



where 



^ := (i + 2- +fc+i r - 1 < I l^Z/ || J < 1; ■ (26) 

In (p5|) , the term involving becomes negligible compared to the other one as 
n — ► oo and we are left with : 

C 2-(«- 1 /2) _ 1 - 2- n (9-V2) 
|/| - 2~ n 2 1 - 2-(9~V2) ' 

The estimate from below is obtained similarly. | 



Remark 4.2 To calculate M q (a(-),l) for a random function a(x) notice, that our 
estimates depend on Nk linearly. Therefore inequalities for the mathematical expec- 
tation of Nk are enough to prove the statement of Lemma \4-l\ for the mathematical 
expectation of M q (a(-) , I) . 

4.2 Estimate of N k . 

Let £u,(x) be a realization of the standard Brownian motion on some finite interval /, 
and let C w be the convex hull of the realization w(y) := + x) dx. We denote 

by S the set of end points of straight segments (corresponding to shocks) of C w . This 
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set is a closed set of zero Lebesgue measure. For any fixed integer k > 0, we introduce 
the following notation : 

N k := #{A, : 2- k - 1 < |A,| < 2- k }, (27) 

where UjAj = I — S and Aj are intervals, whose end points lies in the set S. 

Let us fix two points xi,X2 G I and a small number < i5 < 1. Consider two 
intervals I\ := [x\ — \x% — X\\5, X\ + \x2~ X\\S) and I2 '■— [x^ — \x% — X\\S, X2 + \x2~ Xi\S). 
We shall denote by ((h, h) the indicator function of the event that C w has a straight 
segment whose endpoints lie inside the intervals h,l2, respectively. 

Lemma 4.3 (Sinai) There is a constant Cq = Cq(5) > such that 

Co 1 ^ - X^ 2 < (C(/l, h)) < C,\X 2 - X X \ X ' 2 . 



Lemma 4.4 There is a constant C > such that C 1 2 k/2 < (N k ) < C2 k/2 for all k 
large enough. 

Proof. Fix k and choose a sufficiently large m. Decompose the segment / onto 
equal intervals Ij := [cj, c^), indexed from the left to the right, of length 2~ h /m and 
consider the pairs Jj, Ij such that m/2 < j — i < m + 2. Then 

2- k - 1 < (j - i)2- k /m < 2- k - 2 • 2- k /m, (28) 

which means that for any pair of points x G U and y G Ij we have 2 -fc ~ 1 < \x — y\ < 
2~ k . Therefore, 

N k = Y,C(h,Ij)- (29) 

u 

By Lemma [4 . 3| the mathematical expectation of ( can be estimated as follows 

CV 1 y/l/2 2- k ' 2 < (C(I h Ij)) < C y/l - 2/m 2~ k / 2 . (30) 
Thus, summing over all pairs of indices, we have : 

c -l 2 k/2 < ^ < ^ m 2 2 k 2 -k/2 = c 2 fe/2 ) ( 31 ) 

which completes the proof. | 
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M q as a function of I 



Figure 2: Moments of order q, as labeled, of Lagrangian increment vs. separation 
I for Brownian initial velocity (h = 1/2) at t = 0.5. The simulations used periodic 
conditions of unit size with a mesh e = 2 -22 . Notice the conspicuous, but spurious, 
scaling with unity exponent at small separations. 



Local scaling exponent for h = 1/2 



Figure 3: Local scaling exponent obtained as the logarithmic derivative of M q (l) with 
respect to I. Same conditions as Fig. |2|. 

5 Simulations and spurious scaling regime 

The numerical strategy for solving Burgers equation with fractional Brownian initial 
velocity has been described in detail in Ref. || (Section 5). Let us just recall some 
key points here. In the simulations, it is necessary to introduce both a large-scale and 
a small-scale cutoff. We find it convenient to work with periodic velocity fields and 
to set the spatial period to unity. Then, the scales accessible are clearly restricted to 
the range e < I < 1, where e = 1/N is the inverse of the number of grid points in the 
simulation. In Ref. [[| up to 2 20 grid points were used. For reasons which will become 
clear, we had here to work with even higher resolution, using N = 2 22 . Moments were 
calculated by averaging over space and over about 1200 realizations. The constant C 
appearing in (0) is always taken unity. 

Fig. [| gives log-log plots of the moments M q (l) for h = 1/2 and 17 values of the 
exponent q varying linearly between zero and two. The graphs appear to all have the 
same unity slope at small Z's and a different g-dependent slope at large Vs. Fig. |3] 
shows the logarithmic derivative of M q (l) with respect to I, a measure of the local 
scaling exponent. It is seen that below a separation I of about 10 -4 all the exponents 
r q go to unity, while at very large separations they appear to approach the value r q = q 
(at least for g's up to one). The latter result is an immediate consequence of the fact 
that, at separations much larger than the coalescence length l c , the Lagrangian map 
is very close to the identity. Hence, M q (l) ~ l q . 

We shall now show that, when < q < = h, the former result T q = 1 is 
a numerical artifact due to a spurious discretization effect, affecting all separations 
such that I < l sp ~ te h , where e is the numerical mesh. By (19), for a given I, the 
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Local scaling exponent for h = 3/4 



Figure 4: Same as Fig. |3|, but for h = 3/4. Notice that the region of spurious scaling 
T q = 1 has shifted to smaller separations and that a region with r q = q/h for small g 
becomes visible. 

dominant contribution to M q (l) should come from those shocks with a Lagrangian 
length Aa = 2~ n such that, 

/ ~tT nh ~t(Aa) h . (32) 

Clearly, no shocks can be represented which have Aa < e, where e is the numerical 
mesh. This gives a cutoff in I at 

k P ~te h . (33) 

For / < l sp , the numerically measured Lagrangian increment will typically take only 
two values: e with probability l/l sp and zero with probability 1 — l/l sp . Hence, 

M q {l) ~eH/l spj (34) 

which implies T q = 1. This is the spurious scaling announced. 

Thus, the scaling regime with r q = q/h, discussed in the theoretical part of this 
paper, can be observed only at scales such that 

Z Bp ~te fc <Z«:Z c ~t 1 /( 1 - h >, (35) 

where 

l c ~t 1/( ±- h \ (36) 

is the coalescence length ||. 

In practice fl3~5j) is a strong constraint : if we also want to avoid contaminations of 
T q due to finite box size, we should take l c ~ significantly smaller than unity. 

Hence, for h = 1/2, the range of nonspurious separations defined by (p5|) will be too 
small to be visible unless we work at extremely high resolutions (very small e's). 

Inspection of ( p5|) reveals, however, another strategy : we can increase h and 
thereby push the spurious range to smaller separations. Since the fractional Brownian 
properties of the velocity disappear at h — 1, we chose h = 3/4 as a trade-off. 
Fig. f| shows the same result as in Fig. |3] but, now, for h = 3/4. We observe that 
the region of spurious scaling has now been pushed below / ~ 10 -5 and that for 
0<g<g* = /i = 3/4a kind of plateau near r q = q/h is seen at the smallest 
separations not affected by spurious scaling. Measuring directly the exponents r q in 
this range produces Fig. |5| where r q is plotted vs. the exponent q, both for h = 1/2 
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Phase transition 



Figure 5: Scaling exponent measured at the smallest separation not affected by spu- 
rious scaling, plotted vs. the order q for h = 1/2 and h = 3/4, as labeled. The 
piecewise linear graphs represent theoretical predictions. 

and h — 3/4. The comparison with our theoretical prediction of Section |3] (the thick 
straight lines) is now satisfactory, the only remaining discrepancies being caused by 
unavoidable finite size effects which soften the phase transition. 
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